Measles Example
Note: This tutorial assumes that you have successfully installed ABSEIR, and are at least passingly familiar with compartmental models. Some introductory information is available in this vignette.
Step 0: Setup
Data for this example were obtained from the measlesWeserEms data set from the surveillance package.
# load the ABSEIR library
library(ABSEIR)## Loading required package: Rcpp## Loading required package: parallel## Loading required package: compilerlibrary(surveillance)## Loading required package: sp## Loading required package: xtable## Loading required package: polyCub## This is surveillance 1.13.0. For overview type 'help(surveillance)'.# Read in data
data(measlesWeserEms)
# Identify cases
cases<-measlesWeserEms@observed
epidemic.start = min(which(apply(cases, 1, max) > 0))
cases = cases[(epidemic.start-1):nrow(cases),]
# Obtain distance matrix
neighbourhood<-measlesWeserEms@neighbourhood
# Week
week <- measlesWeserEms@epoch[(epidemic.start-1):length(measlesWeserEms@epoch)]
# Vaccination and population data set
vaccine.data <- measlesWeserEms@map@data
vaccine.data$adminID <- rownames(vaccine.data)
# Population size
N <- vaccine.data$POPULATION
# Check that spatial unit ordering makes sense
if (!all(vaccine.data$adminID == colnames(neighbourhood))){
stop("Error, make sure spatial unit ordering is consistent.")
}Step 1: Data Model and Initial Values
# Make cumulative
weser.data_model = DataModel(Y=apply(cases, 2,cumsum),
type = "identity", # Assume data is correct
compartment = "I_star", # Data related to new infections
cumulative = TRUE # Not reported on cumulative scale
)
I0 = (apply(cases[1:3,], 2, max) > 0)*2
E0 = I0
R0 = 0*I0
S0 = N-E0-I0-R0
weser.initial_values = InitialValueContainer(S0 = S0,
E0 = E0,
I0 = I0,
R0 = R0)
# No reinfection
weser.reinfection_model = ReinfectionModel("SEIR")Step 2: Exposure Process
To specify the exposure process, we need to create the exposure design matrix. This matrix of temporal and spatial covariates must have a particular structure. Each of the \(J\) spatial locations receives a \(T\times p\) design matrix, where \(p\) is the number of covariates. These matrices are then ‘stacked’ row-wise in the same spatial ordering as presented in the columns and rows of the neighborhood matrix (defined above).
Here, we include the following covariates:
- An intercept term
- The population density
- The proportion of students with a vaccine card
- The proportion of students receiving at least one vaccination
- The proportion of students who were doubly vaccinated
- A three degree of freedom trigonometric temporal basis, to capture seasonality
For more information on these covariates, see the documentation for the ‘measlesWeserEMS’ data.
As you can see, while ABSEIR doesn’t directly fit SEIRV models (which incorporate a vaccination compartment), we can include vaccination data into the exposure probability structure, effectively giving vaccinated persons a different probability of contracting the pathogen of interest.
n.locations <- ncol(neighbourhood)
n.timepoints <- length(week)
time_invariant_covariates <- data.frame(intercept = rep(1, nrow(vaccine.data)),
popDensity = vaccine.data$POPULATION/vaccine.data$AREA,
proportionVaccineCard = vaccine.data$vaccdoc.2004,
proportionVaccine1 = vaccine.data$vacc1.2004,
proportionVaccine2 = vaccine.data$vacc2.2004)
time_varying_covariates <- data.frame(sin_component = sin(week/52*2*pi),
cos_component = cos(week/52*2*pi),
trig_interact = sin(week/52*2*pi)*cos(week/52*2*pi))
exposure.design.matrix <- as.matrix(
cbind(
time_invariant_covariates[rep(1:n.locations, each = n.timepoints),],
time_varying_covariates[rep(1:n.timepoints, n.locations),]
)
)
## Build the exposure model
weser.exposure_model <- ExposureModel(X = exposure.design.matrix,
nTpt = n.timepoints,
nLoc = n.locations,
betaPriorPrecision = rep(1, ncol(exposure.design.matrix)),
betaPriorMean = rep(0, ncol(exposure.design.matrix)))## [1] "Assuming equally spaced count data."Step 3: Contact Structure
We now need to specify the contact structures. Here, we’ll build both a gravity model (which depends on population size and distance), and a simple CAR model (which just depends on shared borders between regions).
# Build a gravity model, in which the contact process between a pair of
# spatial locations is proportional to the product of their populations divided
# by the squared 'distance'
pop.matrix <- matrix(N, nrow = length(N), ncol = length(N))
gravityModel <- (pop.matrix * t(pop.matrix))/neighbourhood^2
diag(gravityModel) <- 1
# Rescale
maxRowSum <- max(apply(gravityModel,1,sum))
gravityModel <- gravityModel/maxRowSum
weser.distance_model <- DistanceModel(list(gravityModel),
priorAlpha = 1,
priorBeta = 1)
# Build a simpler contact model, in which the contact probability between
# a pair of spatial locations is only nonzero when the distance is equal to 1
# (CAR specification)
weser.CAR_model <- DistanceModel(list((neighbourhood == 1)*1),
priorAlpha = 1,
priorBeta = 1
)Step 4: Specify Latent and Infectious Period Transition Models
# 9-12 day latent period
# Infectious
weser.transition_priors = ExponentialTransitionPriors(p_ei = 0.8, # Guess at E to I transition probability (per week)
p_ir = 0.8, # Guess at I to R transition probability (per week)
p_ei_ess = 10, # confidence
p_ir_ess = 10 # confidence
)Step 5: Configure ABC Sampler
weser.sampling_control <- SamplingControl(seed=123124,
n_cores = 14,
algorithm = "Beaumont2009",
params = list(
batch_size = 2000,
init_batch_size = 1000000,
epochs = 1e6,
shrinkage = 0.99,
max_batches = 200,
multivariate_perturbation=FALSE
))Step 6: Run Our Two Models
# Cache the model objects so we don't need to re-run the analysis for things like format tweaks
if (!file.exists( "weserModel1.rda")){
t1 <- system.time(
{
weser.model1 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.distance_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control,
samples = 50,
verbose = FALSE)
}
)
save("weser.model1", "t1", file = "weserModel1.rda")
} else {
load("weserModel1.rda")
}
if (!file.exists( "weserModel2.rda")){
t2 <- system.time(
{
weser.model2 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.CAR_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control,
samples = 50,
verbose = 2)
}
)
save("weser.model2", "t2", file = "weserModel2.rda")
} else {
load("weserModel2.rda")
}
print(t1)## user system elapsed
## 10504.509 16.417 760.814print(t2)## user system elapsed
## 23474.217 34.397 1701.912Parameter Summaries
summary(weser.model1)## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -0.350 0.719 -1.752 0.946
## Beta_SE_2 0.475 0.790 -1.299 1.995
## Beta_SE_3 -0.632 0.615 -1.944 0.497
## Beta_SE_4 0.108 0.518 -0.966 0.923
## Beta_SE_5 0.708 0.890 -1.123 2.057
## Beta_SE_6 0.050 0.203 -0.359 0.414
## Beta_SE_7 0.224 0.318 -0.466 0.830
## Beta_SE_8 0.124 1.570 -3.647 2.266
## rho_1 0.116 0.036 0.043 0.183
## gamma_EI 1.275 0.660 0.244 2.341
## gamma_IR 1.574 0.534 0.757 2.575summary(weser.model2)## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -0.201 0.307 -0.792 0.495
## Beta_SE_2 0.103 1.181 -2.181 2.461
## Beta_SE_3 -0.176 0.386 -0.980 0.509
## Beta_SE_4 -1.018 0.578 -1.833 0.084
## Beta_SE_5 0.921 0.331 0.300 1.708
## Beta_SE_6 0.162 0.425 -0.648 0.883
## Beta_SE_7 -0.378 0.406 -1.008 0.349
## Beta_SE_8 -0.139 0.868 -1.772 1.392
## rho_1 0.026 0.010 0.009 0.044
## gamma_EI 1.595 0.396 0.728 2.179
## gamma_IR 1.280 0.754 0.417 2.709Graphical Summary Function
This function will plot cases each of the spatial locations for a model, add a black line for the mean posterior predictive number of cases, and add blue dashed lines for the 99% credible interval for the number of cases.
makePlots <- function(modelObj, nm){
sims <- epidemic.simulations(modelObj, replicates = 50)
Is <- lapply(sims$simulationResults, function(x){x$I_star})
Is <- array(Reduce(c, Is), dim = c(nrow(Is[[1]]),
ncol(Is[[2]]),
length(Is)))
Ism <- apply(Is, 1:2, mean)
Islb <- apply(Is, 1:2, quantile, probs = c(0.025))
Isub <- apply(Is, 1:2, quantile, probs = c(0.975))
plotLocation <- function(x, model){
plot(cases[,x], ylim = c(0, max(Isub[,x])),
main = paste(model, ": Observed and Posterior Predictive Simulation.\n location ",
colnames(neighbourhood)[x], sep = ""))
lines(Ism[,x], col = rgb(0,0,0,0.8), lwd = 2)
lines(Islb[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
lines(Isub[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
#apply(Is, 3, function(i){
# lines(i[,x], col = rgb(0,0,0,0.1))
#})
points(cases[,x], pch = 16, col = "blue")
}
for (i in 1:ncol(neighbourhood)){
plotLocation(i, nm)
}
}Graphical Illustration of Model 1:
makePlots(weser.model1, "Dist")
Graphical Illustration of Model 2:
makePlots(weser.model2, "CAR")
Bayes Factors Comparing the Models (row v. column):
comps <- compareModels(modelList = list(weser.model1, weser.model2), n_samples = 100)## [1] "Assuming equal prior probabilities."rownames(comps) <- colnames(comps) <- c("Distance", "CAR")
print(comps)## Distance CAR
## Distance 1.0 0.1315789
## CAR 7.6 1.0000000It looks like the preferred model is the CAR model.